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We present a detailed description of semi-quantum molecular dynamics simulation of stochastic 
dynamics of a system of interacting particles. Within this approach, the dynamics of the system is 
described with the use of classical Newtonian equations of motion in which the effects of phonon 
quantum statistics are introduced through random Langevin-like forces with a specific power spec- 
tral density (the color noise). The color noise describes the interaction of the molecular system with 
the thermostat. We apply this technique to the simulation of thermal properties and heat trans- 
port in different low-dimensional nanostructures. We describe the determination of temperature in 
quantum lattice systems, to which the equipartition limit is not applied. We show that one can 
determine the temperature of such system from the measured power spectrum and temperature- and 
relaxation-rate-independent density of vibrational (phonon) states. We simulate the specific heat 
and heat transport in carbon nanotubes, as well as the heat transport in molecular nanoribbons 
with perfect (atomically smooth) and rough (porous) edges, and in nanoribbons with strongly an- 
harmonic periodic interatomic potentials. We show that the effects of quantum statistics of phonons 
are essential for the carbon nanotube in the whole temperature range T < 500K, in which the val- 
ues of the specific heat and thermal conductivity of the nanotube are considerably less than that 
obtained within the description based on classical statistics of phonons. This conclusion is also ap- 
plicable to other carbon-based materials and systems with high Debye temperature like graphene, 
graphene nanoribbons, fuUerene, diamond, diamond nanowires etc. We show that the existence 
of rough edges and quantum statistics of phonons change drastically the low temperature thermal 
conductivity of the nanoribbon in comparison with that of the nanoribbon with perfect edges and 
classical phonon dynamics and statistics. The semi-quantum molecular dynamics approach allows 
us to model the transition in the rough-edge nanoribbons from the thermal-insulator-like behavior 
at high temperature, when the thermal conductivity decreases with the conductor length, to the 
ballistic-conductor-like behavior at low temperature, when the thermal conductivity increases with 
the conductor length. We also show how the combination of strong nonlinearity of periodic inter- 
atomic potentials with the quantum statistics of phonons changes completely the low-temperature 
thermal conductivity of the system. 



I. INTRODUCTION 

Molecular dynamics (MD) is a method of numerical 
modeling of molecular systems based on classical Newto- 
nian mechanics. It does not allow for the description 
of pure quantum effects such as freezing out of high- 
frequency oscillations at low temperatures and the re- 
lated decrease to zero of heat capacity for T ^ 0. In 
classical molecular dynamics (CMD), each dynamical de- 
gree of freedom possesses the same energy ksT, where ks 
is Boltzmann constant. Therefore, in classical statistics 
the specific heat of a solid almost does not depend on 
temperature when only relatively small changes, caused 
the anharmonicity of the potential, can be taken into 
account [1]. On the other hand, because of its complex- 
ity, a pure quantum-mechanical description docs not al- 
low in general the detailed modeling of the dynamics of 
many-body systems. To overcome these obstacles, dif- 
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ferent semiclassical methods, which allow to take into an 
account quantum effects in the dynamics of molecular 
systems, have been proposed [2-8]. The most convenient 
for the numerical modeling is the use of the Langevin 
equations with color- noise random forces [5, 7]. In this 
approximation, the dynamics of the system is described 
with the use of classical Newtonian equations of motion, 
while the quantum effects arc introduced through ran- 
dom Langevin-like forces with a specific power spectral 
density (the color noise), which describe the interaction 
of the molecular system with the thermostat. Below we 
give a detailed description of this semi-quantum molec- 
ular dynamics (SQMD) approach in application to the 
simulation of specific heat and heat transport in differ- 
ent low-dimensional nanostructures. 

The paper is organized as follows. In Section II we 
describe the temperature-dependent Langevin dynam- 
ics of the system under the action of random forces. 
If the random forces are delta-correlated in time do- 
main, this corresponds to the white noise with a flat 
power spectral density. This situation corresponds to 
high enough temperatures, when /cbT is larger that 
the quantum of the highest-frequency mode in the sys- 
tem, ksT ^ hujm. But for low enough temperature. 



2 



fcsT <C huJm, the stochastic dynamics of the system 
should be described with the use of random Langevin- 
Uke forces with a non-flat power spectral density, which 
corresponds to the system with color noise. In Section 
III we describe the determination of temperature in quan- 
tum lattice systems, to which the equipartition limit is 
not applied. Wc show how one can determine the tem- 
perature of such system from the measured power spec- 
trum and temperature- and relaxation-rate-independent 
density of vibrational (phonon) states. In Section IV we 
describe a method for the generation of color noise with 
the power spectrum, which is consistent with the quan- 
tum fluctuation-dissipation theorem [1, 9]. In Sections 
V and VII we apply such semi-quantum molecular dy- 
namics approach to the simulation of the specific heat 
and phonon heat transport in carbon nanotubcs. In Sec- 
tion VI we apply this method to the simulation of heat 
transport in a molecular nanoribbon with rough edges. 
Previously, we have predicted analytically [10] and con- 
firmed by classical molecular dynamics simulations [11] 
that rough edges of molecular nanoribbon (or nanowire) 
cause strong suppression of phonon thermal conductiv- 
ity due to strong momentum-nonconserving phonon scat- 
tering. Here we show that the rough edges and quan- 
tum statistics of phonons change drastically the thermal 
conductivity of the rough-edge nanoribbon in compari- 
son with that of the nanoribbon with perfect (atomically 
smooth) edges. The semi-quantum molecular dynamics 
approach allows us to model the transition in the rough- 
edge nanoribbons from the thcrmal-insulator-likc behav- 
ior at high temperature, when the thermal conductivity 
decreases with the conductor length, see Ref. [11], to 
the ballistic-conductor-likc behavior at low temperature, 
when the thermal conductivity increases with the con- 
ductor length. In Section VIII we apply this technique 
to the simulation of thermal transport in a nanoribbon 
with strongly anharmonic periodic interatomic poten- 
tials. Wc show how the combination of strong nonlin- 
earity of the interatomic potentials with quantum statis- 
tics of phonons changes completely the low-tcmpcrature 
thermal conductivity of the nanoribbon. In Section IX 
we provide a summary and discussions of all the main 
results of the paper. 



II. LANGEVIN EQUATIONS WITH COLOR 
NOISE 



In the presence of a Langevin thermostat, equations of 
motion of coupled particles have the following form: 



d 

M„r„ = -T^H - M„rr„ + S„ 



(1) 



where the three-dimensional radius- vector r„(i) gives 
the position of the n-th particle at the time instant t 
{n = 1,2, . . . , N, N being the total number of particles), 
M„ is the particle mass, H is the Hamiltonian of the sys- 
tem and F„ = —dH/dVn gives the force applied to the 



n-th particle caused by the interaction with the other 
particles, F = 1/t,. is the friction coefficient (<r being the 
relaxation time due to the interaction with the thermo- 
stat), and S„ = {Cct, )i}a=i are random forces with the 
Gaussian distribution. 

In the semi-quantum approach the random forces do 
not represent in general white noise. The power spectral 
density of the random forces in that description should 
be given by the quantum fluctuation-dissipation theorem 
[1, 9]: 



2MnTkBTdal3Sn7nP{'^, T) , 



(2) 



where n,m = 1, 2, N , a, P = 1, 2, 3 are Cartesian com- 
ponents, 5nm and Safi are Kronecker delta-symbols. Here 
the positive and even in frequency oj function 
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gives the dimensionless power spectral density at temper- 
ature T of an oscillator with frequency lo. The first term 
in Eq. (3) gives the contribution of the zero-point oscil- 
lations to the power spectral density of random forces. 
With the use of Eqs. (2) and (3), we get the correlation 
function of the random forces as 
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(^anWC/3m(0)) = 

p{u;,T)e'-^*^ . (4) 



Let flmax be the maximal vibrational eigenfrequency 
of the molecular system. At high temperatures, T ^ 
ft^max/kB, one has p{uj,T) = 1 and the molecular sys- 
tem will perceive the random forces as a delta-correlated 
white noise: 



{(,an{t)£,l3m{0)) = 2AInT kBTSal3SnviS{t) . 



(5) 



Therefore for high enough temperatures we deal with 
molecular dynamics with Langevin uncorrelated random 
forces. But for low temperatures, the random forces are 
correlated according to the function given by Eqs. (3) 
and (4). In other words, for low temperatures the ran- 
dom forces in the system represent the color noise with 
the dimensionless power spectral density p{uj, T) given by 
Eq. (3). 



III. DETERMINATION OF TEMPERATURE IN 
QUANTUM LATTICE SYSTEMS 

In connection with the consideration of the thermo- 
dynamics and dynamics of non-classical lattices, in this 
Section we discuss possible determination of the temper- 
ature in the systems, to which the equipartition limit is 
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not applied. Wc start from Eqs. (1), rewritten as the 
equations of motion for the lattice vibrational mode with 
eigenfrequency f2„: 

M„Q„ + M„rQ„ + M^nlQn = S„, (6) 

where S„ = {^Q.n}a=i are random forces with the spec- 
tral density given by Eqs. (2) and (3). 
From Eqs. (2), (3) and (6) we obtain 



M„(Q^W = -coth(^) — 



2r 



2p2 



, (7) 



and power spectrum Et{uj) of lattice vibrations with 
3A^ — 6 nonzero eigenmodes: 
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(8) 



We can use {Q^)ui,t and the power spectrum Et{i^), 
Eqs. (7) and (8), to compute the average energy of the 
system. Taking into account that the above quantities 
are even functions of frequency w, for positive a; the av- 
erage energy of the system can be computed as: 



EiT) 



Et{uj) 
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Et{u}) — . 
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(9) 



In the weakly dissipative limit, with F ^ Qn for all f2„, 
which is realized in all the considered system, for positive 
Lj and n„ one has, in the sense of generalized functions: 
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(10) 



(fi2 - u;2)2 + ^2r2 

We can introduce in this limit the density of vibrational 
(phonon) states (for positive frequencies) as: 

3N <,T^^9 3JV 



2Fr!2 



- {ni - uj^y + w2F2 



nY,S{u^-nn). (11) 



n=7 



Then the power spectrum (for positive uj) consists of an 
array of the delta-functions at the system eigenfrequen- 
cies, weighted by temperature, 

Et{u:) = —coth{^j-^)D{u;) 



>a;coth(;^)E'^(^-^^«), (12) 



and the average energy of the system (9) has the following 
form: 



3N 



E{T) = J2 



n=7 



exjp{hQn / ksT) — 1 



(13) 



Here the first term in the sum gives the temperature- 
independent contribution of zero-point oscillations to the 
energy of the system. As one can see from Eqs. (11)-(13), 
in the weakly dissipative limit the density of vibrational 
states, power spectrum and average energy of the system 
do not depend explicitly on the relaxation rate F. 

From Eq. (12) we get the "quantum definition" of 
temperature T through the power spectrum of lattice vi- 
brations at the given frequency: 



T ~ TiLo I ks In 



1 + ^t(w) 



1 - Ariid) 



where 



Ariuj) = nujD{uj)/{2ET{uj)). 



(14) 



(15) 



In the case when one omits the zero-point contribution 
to the spectrum of the random forces, namely considers 



cxjp{huj/kBT) — 1 



(16) 



in Eqs. (2) and (6), the temperature can be determined 
from the following equation: 



T = fiuj/{kB ln[l + 2AT(w)]), 



(17) 



where in the definition of the function At{uj) in Eq. (15) 
the function Et{uj), in contrast to function ET{io) in Eq. 
(8), does not have now the zero-point contribution, 



Et{uj) 



3N 



exp(?ia;/fcBr) 
and, correspondingly, one has 



2T^l 



- t^2)2 + ^2^2 ' 

(18) 
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exp^hfln/kBT) 



1 



(19) 



In the equipartition limit, which is realized for T ^ 
hilmax/kB; one has Et{'-^) = kBTD{uj), At(i^) «C 1, and 
both Eqs. (14) and (17) turn into the identity T = T. 

Therefore one can determine the temperature of 
the quantum lattice system from the measured power 
spectrum Et{ijj) and temperature- and relaxation-rate- 
independent density of vibrational (phonon) states D{lo). 
In Fig. 6 we show the power spectrum, computed within 
the semi-quantum approach with F = l/tr, = 0.4 ps, 
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line 1, and in the limit, when the spectrum is given by 
an array of the smoothed delta-functions, line 3. As one 
can see in this figure, the particular choice of the (long 
enough) relaxation time tr = 1/T indeed does not affect 
the presented results. The determination of tempera- 
ture, given by Eqs. (14) and (17), is valid for weakly 
anharmonic systems, in which the power spectrum of 
lattice vibrations is close to the one, given by Eq. (8). 
Carbon-based materials like carbon nanotubes, graphene 
and graphene nanoribbons, and crystal structures with 
stiff valence bonds belong to such systems. It is worth 
mentioning that the temperature, given by Eq. (14) or 
Eq. (17), does not depend on frequency uj only in the 
case of equilibrium harmonic lattice excitations, driven 
by random forces with power spectrum given by the Bose- 
Einstein distribution, Eqs. (2) and (3). The anharmonic- 
ity of the lattice can change the power spectrum of the ex- 
citations, see Fig. 12 (b) below. In such case the "phonon 
temperature" starts to depend on frequency (similar to 
the frequency- and direction-dependent temperature of 
photons in non-equilibrium radiation, see Ref. [1]), and 
"hot phonons" with non-equilibrium high frequencies 
appear in the system, see Section VI below. 



IV. COLOR NOISE GENERATION 



forces in Eqs. (1) as 

= 2M„rfcsT(5„„5^„,)^, (21) 

such that 

f duj 
2^ 
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-M^T{kBTf / {S^nSf„n)^e'^^7^ , (22) 



and 



(23) 



Here, the uncorrelated random functions Son(''") ^^nd 
Sin(T) will generate, in a finite frequency interval, the 
power spectra po{uj) and pi(w), respectively. 

For a molecular system with maximal vibrational 
eigenfrequency Wm, we will construct the dimensionless 
random vector function Son(T), whose power spectral 
density generates Po{u}) in the frequency interval [0, Wm] 
(wm = haJm/ksT). Each component of this vector can 
be written as a sum of dimensionless random functions 



For the implementation of the semi-quantum approach 
in molecular dynamics, random forces with the power 
spectral density given by p{uj) must be generated. Ex- 
isting numerical techniques allow the generation of a 
random quantity from a prescribed correlation function 
[12]. This approximation was used in [7] in the modeling 
of thermodynamic properties of liquid ^He above the A 
point. But universal methods require in general the usage 
of high numerical resources. For instance, the utilization 
of the method proposed in [12] requires the execution of 
the Fast Fourier Transform for every random force value. 
We propose a technique which takes into account the 
specific form of the dimensionless power spectral density 
p(a;), where uj = fiu /ksT is a dimensionless frequency. 

In terms of d), the dimensionless power spectral density 
of random forces looks as 



piuj) = -w + -— - 

2 exp(w) 



-wcoth(<I>/2). (20) 



The random force with dimensionless power spectral 
density, given by Eq. (20), can be presented as a sum 
of two independent functions, p(oj) = P(){oj) + Pi{'^), 
where the first function, Pq{u)) = w/2, has a linear 
power spectral density, while the second one is pi{u)) ~ 
a'/[exp(w) — 1]. 

To generate the color noise with a given dimensionless 
power spectral density p{u)) , we will construct the dimen- 
sionless random vector functions S„(r) = {Sa^n\^a=i ~ 
Sori(''") + Si„(t) of the dimensionless time r = tksT /Ti, 
which will give the power spectral density of the random 



Soan{T) = ^Ci[77Q„.i(r) - C,ans{T)]. 



(24) 



The functions {CaniliLi which generate the color noise, 
are the solutions of the linear equations: 



Cni{T) = A,[r/„m(T) - Cam(T)], 



(25) 



where Canii'^) derivative with respect to r. Here 

VaniiT) are dimensionless white- noise random forces, 
with the correlation functions 



{llcmi{T)rii3kj{0)) = 2SapSnkStjS{T)/\i 



(26) 



By solving Eqs. (25) in the frequency domain, (ani 
can be substituted in Eq. (24) and the power spectra of 
Son can be obtained as 



(27) 



The dimensionless parameters Ci and can be found by 
minimizing the integral 
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2cfx 



2 ^A.(A? + x2) 



dx 



(28) 



with respect to the parameters ci, Ai, C4, A4, where 
Ci = Ci/uj„i, \i ~ Xi/ujn. The coefficients Cj, Aj, 
i = 1,...,4, obtained within this procedure, are shown 
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TABLE I: Value of the coefficients Ai, Ci, Qt, Fi. 



coefficient 


value 


coefficient 


value 




1.763817 


C5 


1.8315 




0.394613 


C6 


0.3429 


\ / — 
■\3/iOm 


0.103506 


Hs 


2.7189 




0.015873 


He 


1.2223 


Cl/uJm 


1.043576 


fs 


5.0142 


C2/l^m 


0.177222 


fe 


3.2974 




0.050319 








0.010241 
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FIG. 1: (Color online) Power spectral density of color noise 
Sq{t) (blue line). Red line is given by Eq. (27), dashed 
line shows the linear function po(uj)/uJm = uj /2uJrn- The two 
lines practically coincide in the frequency interval 0.015 < 
tj/tJm < 1.1. The inset shows spectral densities given by 
Eq. (27) (red line) and by linear function po(a;)/cjm (dashed 
line) in the frequency interval < uj/uJm < 0.04. For color 
noise generation, Eqs. (25) were numerically integrated with 
the use of fourth order Runge-Kutta method with a constant 
integration step Ar = 0.05/a)m. 



in Tabic I. For these parameters, integral (28) reaches its 
lower value, equal to 2.46 x 10~^. 

As one can sec from Fig. 1, the function given by Eq. 
(27) approximates with high accuracy the linear function 
Po{oj) ~ in the frequency interval 0.015 < uj/uJm < 
1.1. 

On the other hand, the random function 5*10,1 (r), 
which will generate the power spectral density pi(a)), can 
be approximated by a sum of two random functions with 
relatively narrow frequency spectra: 



(29) 



In this sum the dimensionless random functions C,ani {t) , 
1 = 5,6, satisfy the equations of motion as 

Cm(^) = Vani{T) - (^fCamir) - TiCnii'^), (30) 

where, as before, rjaniiT) are delta-correlated white noise 




FIG. 2: (Color online) Power spectral density of color noise 
Si{t), Eq. (32) (blue line). Red line shows the function 
Pi{io) = Lj/[exp{Lj) — 1]. For color noise generation, Eqs. 
(30) were numerically integrated with the use of fourth order 
Runge-Kutta method with a constant integration step Ar = 
0.02. 



functions: 

{VantiT)rii3kj{0)) ^2fiSai3SnkSijS{T). (31) 

We can solve, as previously, Eq. (30) in the frequency 
domain and insert C^ani into Eq. (29) to obtain the power 
spectrum of Sian- 



2c?f. 



^^(n|-^2)2+^2f2 



(32) 



The dimensionless parameters {q, £7^, Fij^^g can be 
found by minimizing the integral 

6 



exp(a;) 



4 = 5 ^ « 



(^1 -Cj2)2 +^2f2 



dui (33) 



with respect to the parameters C5, F5, Sis, ce, Fg, f^g. The 
obtained parameters are shown in Table I. For these pa- 
rameters, integral (33) reaches its lower value, equal to 
2.3 X 10^*. Since the integral (33) is determined by the 
rapid-descending functions, the infinite integration limit 
[0, 00) can be replaced by a finite one [0, Wm]- Numerical 
integration of the integral (33) shows that its minimal 
value practically does not depend on the upper integra- 
tion limit for Com > 20. 

As one can see from Fig. 2, the function, given by 
Eq. (32), approximates with high accuracy the function 
Pi{uj) = <I'/[exp(a}) — 1] in the frequency interval < < 
10. 

Therefore in the semi-quantum molecular dynamics 
approach one has to solve the Langevin equations (1) 
with random forces 5„ = {<^Qri}a=i, whose power spec- 
tral density is determined as: 



2MnTkBT[{SoanSo^rn)u} + (SlanSl^rn) o] 
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FIG. 3: (Color online) Power spectral density of color noise 
So{t) + Si{r), given by the sum of Eqs. (27) and (32) (blue 
line). Red line shows the function p{io) = (cj/2) coth(a}/2). 
For color noise generation, Eqs. (25) and (30) were numer- 
ically integrated with the use of fourth order Runge-Kutta 
method with a constant integration step At — 0.005. Maxi- 
mal dimensionless frequency LJm ~ 10. 



E 



-,2\2 



(34) 



From Eq. (22) we get the relation between the dimension 
ian{t) and dimensionless San(j) random forces: £,an{t) = 

(anihT/kBT) = kBT^2Mr^r/n[Soc,n{T) + 5la„(r)]. 

In Fig. 3 we show the comparison of the power spec- 
tral density for the dimensionless frequency w, given 
by the sum of Eqs. (27) and (32), with the function 
(a;/2) coth(a;/2), Eq. (20). We see a very good coinci- 
dence in the frequency interval [0, ujm] for uJm = 10. 

To describe the dynamics of molecular system with an 
account for quantum statistics of molecular vibrations 
but without an account for zero-point oscillations, one 
has to keep only the last two terms, with i = 5,6, in 
Eq. (34) and to solve numerically Eqs. (30) in order to 
simulate random forces in Eq. (1). Equations (30) play 
the role of the filter, which filters out the dimensionless 
high-frequency (" non-quantum" ) component of the white 
noise at low temperature. 

It is worth mentioning in this connection that the con- 
sidered semi-quantum approach permits one to obtain 
the correct value of the energy of zero-point oscillations. 
A semi-quantum account for zero-point oscillations in the 
modeling of the properties of liquid '*IIe above the A point 
has allowed to describe correctly the liquid state of the 
helium [7], while the classical molecular dynamics de- 
scription (without an account for zero-point oscillations) 
predicts the solid state of the helium at the same low 
temperature. 



FIG. 4: Form of a single-walled carbon nanotube C300 with 
the zig-zag structure, which consists of A'^ = 300 carbon 
atoms. 



V. TEMPERATURE DEPENDENCE OF 
SPECIFIC HEAT OF CARBON NANOTUBE 

For the illustration of the implementation of the semi- 
quantum approach in molecular dynamics, in this Section 
we find the temperature dependence of specific heat of 
a single-walled carbon nanotube. We consider the nan- 
otube with the (5,5) chirality index (the zig-zag struc- 
ture), which consists of = 300 carbon atoms, see Fig. 
4. In the following we will use the molecular-dynamics 
model of the carbon nanotube, which was discussed in 
details in Rcf. [13]. 

To model nanotube dynamics at temperature T in the 
classical approximation, we will solve the Langevin equa- 
tion (1) with a white noise, where H is a Hamiltonian of 
the nanotube and all M„ = M, M being a mass of a 
carbon atom. The correlation functions of the random 
forces S„, which describe the interaction of an n-th atom 
with the thermostat, are given by Eq. (5). 

We take the initial conditions, which correspond to the 
equilibrium stationary state of the nanotube at T = 0, 
and integrate a system of equations of motion (1) for 
the time t = 20t,-, during which the system will reach 
the equilibrium state at finite temperature. The integra- 
tion beyond this time will allow us to find the average 
thermal energy of the system (H) = E{T). Then the 
specific heat of the molecular system can be found from 
the temperature dependence of the average thermal en- 
ergy C{T) = dE{T)/dT. As one can see in Fig. 5, the 
average thermal energy of the carbon nanotube, placed 
in classical Langevin thermostat, is strictly proportional 
to temperature, E{T) = (3iV - &)kBT « 3NkBT, line 1, 
and, correspondingly, nanotube specific heat almost docs 
not depend on temperature, line 3. This result shows 
both the correctness of the modeling and that carbon 
nanotubes are stiff structures which possess weak anhar- 
monicity of the dynamics of the constituting atoms. One 
can also see in Fig. 5 that the average thermal energy 
per mode e(T) without an account for zero-point oscil- 
lations has the feature that e(T) < fcsT in general and 
e(T) ^ fcfiT for temperature much less than the Debye 
one. For instance, in the carbon nanotube e(T) is almost 
six times less than kBT at room temperature T=300 K. 

To obtain power spectra of the thermal oscillations of 
the atoms in the carbon nanotube, one has to switch 
off the interaction with the thermostat after the estab- 
lishment of the thermal equilibrium in the system. In 
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300 400 

r(K) 

FIG. 5: (Color online) Temperature dependence of the aver- 
age thermal energy per mode e — E/3NkB (a) and specific 
heat per mode c = C/3NkB (b) of a single- walled carbon nan- 
otube Caoo- Lines 1 and 3 show the dependencies, obtained 
with classical Langevin thermostat, red circles give the result 
of numerical modeling. Lines 2 and 4 give the dependen- 
cies, obtained within the quantum description with the use of 
eigenmode spectra without an account for zero-point oscilla- 
tions, Eqs. (19) and (35), blue circles give the result obtained 
within the semi-quantum molecular dynamics approach. 



other words, one has to solve numerieally the frictionless 
equations of motion, Eqs. (1) with F = and S„ = 0, 
with the initial conditions for r„ and r„, which corre- 
spond to the thermalized state of the molecular system. 
By performing the Fast Fourier Transform of the time- 
dependent particle velocities r„ (t) , one can get the power 
spectrum Et{(^), Eq. (8). To increase the accuracy of 
the measurement, it is necessary to perform the averag- 
ing of the obtained results on different initial thermalized 
states of the system. Power spectral density of the ther- 
mal vibrations of the carbon nanotube atoms is shown 
in Fig. 6. As one can see from this figure, all vibrational 
modes of the system are excited in the classical Langevin 
thermostat. 

Carbon nanotube is a rigid structure in which the an- 
harmonicity of atomic dynamics at room temperature is 
weakly pronounced: as one can see from line 2 in Fig. 6, 
the characteristic Debye temperature of the carbon nan- 
otube is very high, ~ huJm/kB ~ 2400K, where 
LOrri ~ 1640 cm~^ is the maximal phonon frequency in 
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FIG. 6: (Color online) Power spectral density of thermal 
vibrations of carbon nanotube C300 atoms at temperature 
T = 300K. Lines 1 (blue) and 2 (red) give the densities ob- 
tained within the semi-quantum and classical descriptions, 
respectively. Line 3 (black) shows the density obtained with 
the use if Eq. (37) for the system of quantum oscillators, 
where the discreteness of the spectrum was smoothed for con- 
venience. The frequency spectrum of the mean kinetic energy 
per atom is shown. 



the system. Therefore we can obtain the temperature de- 
pendence of the carbon nanotube specific heat with the 
use of the spectrum of the harmonic (non-interacting) 
vibrational eigenmodes: C(T) = dE{T)/dT, where the 
average energy of the system E{T) is given by Eq. (13). 
To obtain the eigenfrequencies we have to find all the 
eigenvalues of the square symmetric matrix d^H/dVndrm 
of the 3A^ X 3A^ size: if A„ is one of the eigenvalues, the 
eigenfrequency is determined as r2„ = A„/M„. For 
N = 300 atoms in the nanotube, we have 3iV = 900 
eigenmodes, from which the first 6 modes, which describe 
rigid translations and rotations of the nanotube, have 
zero frequency, ili = . . . = Qq = 0, while the other eigen- 
frequencies are nonzero and are distributed in the interval 



25.1 cm"^ < fir, < 1640.0 cm" 



7, 8,..., 900. We 



can obtain the temperature dependence of the nanotube 
specific heat with the use of Eq. (13): 



3N 



C{T) = kBY. 



( hfln 



^^\kBTj [exp(na./fci3T) - 1]2 



(35) 



As one can see in Fig. 5, the nanotube specific heat in 
the quantum description monotonously goes to zero for 
T ^ 0, and the effects of quantum statistics of phonons 
are essential for the carbon nanotube in the whole tem- 
perature range T < 500K, in which the specific heat C(T) 
of the nanotube is considerably less than the classical- 
statistics value (37V — 6)kB- This conclusion also applies 
to other carbon-based materials and systems with high 
Debye temperature like graphene, graphene nanoribbons, 
fuUerene, diamond, diamond nanowires etc. 

To model the nanotube stochastic dynamics in the 
semi-quantum approach, we will use the Langevin equa- 
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tions of motion (1) with random forces 2ji = {£,cm}a=i 
with the power spectral density, given by {£,anS.i3m)ui = 
2MrkBT{SianSipm)<u, See Eq. (34). Here it is exphc- 
itly taken into account that zero-point osciUations do not 
contribute to the specific heat of the system, and there- 
fore one needs to account only for the random functions, 
given by Eq. (32). 

Therefore the semi-quantum approach in this case 
amounts to the integration of Eqs. (1) and (30), to 
find the average energy (H) = E{T) and the specific 
heat C{T) ~ dE{T)/dT of the system. As one can see 
from Fig. 5, the temperature dependence of the carbon 
nanotube specific heat, obtained by studying the nan- 
otube stochastic dynamics with color noise, coincides 
completely with the dependence, obtained from the spec- 
trum of the harmonic vibrational eigenmodcs in the sys- 
tem. Power spectral density of carbon atoms vibra- 
tions, see Fig. 6, shows that the semi-quantum approach 
describes the thcrmalization of only the low-frequency 
eigenmodes, with oj < lot = ksT/h. In such a way, 
the approach models the quantum freezing of the high- 
frequency eigenmodes at low temperature T < Td. For 
the carbon nanotube, the classical description is defi- 
nitely not valid at room temperature T = 300K, when the 
characteristic frequency wt — 208.5 cm~^ is much less 
than the maximal one in the system, Wm = 1640 cm~^, 
see Fig. 6. 

It is worth mentioning in this connection that the ob- 
tained results for the temperature dependence of the spe- 
cific heat of the carbon nanotube C300 almost do not de- 
pend on the relaxation time t^- All the five values of tr, 
U — 0.1, 0.2, 0.4, 0.8, 1.6 ps, give the same average en- 
ergy E{T) and specific heat C(T) of the nanotube. One 
cannot use either too short tr, when the system dynamics 
becomes the forced one, or too long tr, which increases 
the computation time. We find that the optimal relax- 
ation time for the carbon nanotube is tr = 0.4 ps. 

Above we have introduced and discussed the lattice 
temperature of the system embedded in the Langevin 
thermostat with color noise determined by the quan- 
tum fluctuation-dissipation theorem, see Eqs. (14) and 
(17). But to compute the thermal conductivity, one also 
needs to determine the local temperature of the part of 
the molecular system, which is placed between two ther- 
mostats with different temperatures. In the classical- 
statistics approximation, the lattice temperature can be 
determined as the average double kinetic energy per 
mode (or per degree of freedom): 



T : 



3Nkf 



N 

(E 

n=l 



(36) 



where is a number of atoms and is assumed that 
3A — 6 w 3A^. As we have shown in Chapter III, such 
determination of temperature does not work in the semi- 
quantum approximation when one needs to deal with the 
density and population of vibrational (phonon) states. In 
this case the determination of temperature is given by 



Eq. (14) or by Eq. (17). But in the MD simulations, one 
can also use the determination of temperature, which is 
equivalent to the integral presentation of Eq. (12). 

Having in mind further simulation of thermal conduc- 
tivity, we omit in Eqs. (2) and (6) the zero-point contri- 
bution to the spectrum of the random forces. With the 
use of Eqs. (9), (12), (18) and (19), we find the average 
energy per mode (per degree of freedom) as 



e{T) 



1 

37V 

SAT 



E{T) 



1 

3iV 



Mir 



—y- 

3N ^ exp{hnn/kBT) - 1 ' 



^Dicu) — 

i TT 

(37) 



where D{oj) is density of vibrational states, cf. Eq. (11). 

From the MD simulations, we can also measure the 
spectral density efci„(aj) of the average double kinetic en- 
ergy per mode Ckin as 



N 



dui 



(38) 



In the case of correct " quantum" level occupation, the 
both Eqs. (37) and (38) should correspond to the same 
value of temperature, therefore we obtain the following 
equation for the (numerical) determination of the lattice 
temperature: 



e(T) = eki 



(39) 



Since function e{T) monotonously increases with T, Eq. 
(39) has a unique solution for the temperature. Similar 
equation, but with an account for zero-point oscillations, 
for the rcscaling of MD temperature for an effective one 
was used in Ref. [14] for quantum corrections of MD re- 
sults. As we have mentioned above in discussion of Fig. 5, 
the neglect of zero-point oscillations in the computation 
of the average thermal energy leads to the feature that 
Gfem £ kBT in general and euin ^ for temperature 
much less than the Debye one, see also Fig. 11 (b) and 
(c) below. 

Essentially one can use Eq. (39) for the determination 
of the lattice temperature only for the quantum-statistics 
population of all phonon states. But the situation can 
emerge when high-frequency phonons have an excess ex- 
citation, in comparison with the Bose-Einstein distribu- 
tion for the given temperature, see section VII below. In 
this case the lattice temperature T can be determined 
with the use of only part of the phonon spectrum: 



1 huj 



ekin (w) 



duj 



(40) 

where [cL'i,a;2] is the frequency interval with the Bose- 
Einstein distribution of the average energy per mode 
e(a;). For lui = 0, cj2 > ^max + T and fc^T ^ huj2, 
Eq. (40) gives the classical definition of temperature, 
Skin = kBT. Clearly Eq. (39), which is more convenient 
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for the numerical modeling, gives the correct value of 
lattice temperature only for = and UI2 > ^max + T. 

If the molecular system is completely embedded in 
the Langevin thermostat with color-noise random forces, 
their spectrum, given by Eq. (3), provides the correct 
"quantum" population of all the phonon states. As one 
can see in Fig. 6, at T = 300K the spectral density of the 
average energy per mode efci„(w) in the C300 carbon nan- 
otube (red line 1) coincides exactly with the smoothed 
function D{uj)hu}/{cxp{huj/kBT) - 1)/3N (black line 3). 
In this case Eq. (39) determines the temperature, which 
is exactly equal to the one of the thermostat. 



VI. THERMAL CONDUCTIVITY OF 
NANORIBBON WITH ROUGH EDGES 

First we apply the semi-quantum approach for the 
molecular dynamics simulation of thermal conductivity 
in the nanoribbon with rough edges and harmonic inter- 
particle potential. In such system the vibrational eigen- 
modes do not interact and the considered semi-quantum 
approach turns to be the exact one. 

It was predicted analytically and confirmed by clas- 
sical molecular dynamics simulations that rough edges 
of molecular nanoribbon (or nanowire) cause strong sup- 
pression of phonon thermal conductivity due to strong 
momentum-nonconserving phonon scattering [10, 11]. In 
the case of nonlinear interatomic potential, the ribbon 
has a finite, length-independent thermal conductivity, 
while in the case of harmonic interatomic potential, the 
thermal conductivity decreases with the ribbon length 
(and the ribbon behaves as a thermal insulator). Essen- 
tially for both interatomic potentials, the thermal con- 
ductivity increases with the length of the nanoribbon 
with perfect (atomically smooth) edges. In the case of 
harmonic interatomic potential, phonons in the nanorib- 
bon with rough edges experience effective localization 
due to strong antiresonance multichannel reflection from 
side atomic oscillators in the rough edge layers (the 
Anderson-Fano-like localization due to interference ef- 
fects in phonon backscattcring), see [10, 11]. Apparently 
such strong backscattering can localize phonons with the 
wavelength shorter than the ribbon length. Within the 
classical description, such effect does not depend on tem- 
perature. But this picture will be changed if we take 
into account the quantum effects. The latter result in 
thermalization and correspondingly in participation in 
low temperature thermal transport of long wave phonons 
only. The long wave phonons are much less affected 
by surface roughness than the short wave phonons and 
therefore the effect of thermal conductivity suppression 
should decrease with the decrease of temperature. Be- 
low we demonstrate this effect within the proposed semi- 
quantum approach. 

We consider the system which consists of K parallel 
molecular chains in one plane [11]. Let k be the chain 
number, n be the molecular number in the chain, then the 




FIG. 7: Molecular ribbon made of K — 12 molecular chains. 
Density and widths of rough edges are d — 0.95 and Ki = 4, 
respectively. Lines connect the interacting atoms. 



equilibrium position of the n-th atom in the fc-th chain 
wiU he = na + [I + (— l)'^]a/4, y° „ = bk, where a 
and b are, respectively, the intra- and interchain spatial 
periods, see Fig. 7. 

Only the longitudinal displacements are taken into ac- 
count in a simplest scalar model of two-dimensional crys- 
tal: Xk,n{t) = xl^^+auk,nit), yk,n = yfe_„, wherc Uk^ is a 



dimensionless longitudinal displacement of the (fc, n)-th 
molecule from its equilibrium position. Hamiltonian of 
the system has the form as 



K N 



k=l Tl=l 

K N-1 K 
+ EoY^Y. y{^Kn+l - Uk,n) + So ^ -Bfc, (41) 



fe=l ri=l 



k=l 



where is a number of molecules in each chain, M is 
a mass of molecule, Eq is a characteristic interaction en- 
ergy, V{p) is the dimensionless potential of the nearest- 
neighbor intrachain interactions, p describes relative dis- 
placements of the nearest-neighbor atoms, Ek describes 
the interchain interaction. 

Dimensionless energy of the nearest-neighbor inter- 
chain interactions for odd k is 



N 



Ek ^ U{uk,n - Wfc+l,„-l) -I- 2_^U{uk+l,n - Uk,n), 



n=2 

and for even k is 

N 



N 
n=l 



(42) 



-Efc = ^ U{uk,n - Uk+l,n) + ^ U{uk+l.n+l 



■ Uk,n) 



n=l 



(43) 

where U(p) is the dimensionless potential of the nearest- 
neighbor interchain interaction. 

We will consider a finite rectangular (l<n<A^, 1< 
k < K) with free edges, and harmonic intra- and inter- 
chain interaction potentials: 



V{r) = U{r) = r'^/A. 



(44) 



For the convenience of the modeling, we introduce the 
dimensionless quantities: temperature T ~ T/Tq, time 
i ~ t/to, and energy E = H/Eq, where Tq = It? /Ma?kB, 
to = h/ksTo, and Eq = Ma? /tf^. Then the dimensionless 
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FIG. 8: (Color online) Distribution of (a) local energy flux J, 
(b) local kinetic energy per mode ekin and (c) local tempera- 
ture T along the rough-edge ribbon (ribbon length A'' — 400, 
width K — 12, rough edges widths Ki — 4, porosity of rough 
edges p — 0.05). Gray arias indicate the ribbon ends, em- 
bedded in semi-quantum Langevin thermostats with temper- 
atures T± = 1 ±0.1. 



Hamiltonian will have the form as 



+ 00 



{'>lki(,s)r],nniO)) = '2,TTSkmSin I p(a), T) exp(-iws) ^ , 

J-oo 

(47) 

and dimensionless power spectral density 

p{Q,f) = k^/f)+ ""J^ (48) 
^ exp(a;/ J ) — 1 

where Co = ujIq is the dimensionless frequency. In the 
following we will deal only with the dimensionless quan- 
tities. 

Phonons do not interact in the system with harmonic 
interparticle potential. Therefore the zero-point oscil- 
lations will not affect the thermal transport along the 
nanoribbon and can be not taken into the account. In 
the latter case we can consider the dimensionless ran- 
dom forces with the power spectra, given by {rjknVim) ui = 
2ff{SiknSiim)a, cf. Eq. (34). Random forces in Eq. 

(46) are mn[i) = mnir/f) = fV2rSikn{T). 

Within this approach, we will describe the tempera- 
ture dependence of the nanoribbon thermal conductiv- 
ity k{T). We embed the left end of the ribbon with 
length Ne = 100 in the thermostat with temperature 
r+ = I.IT, while the right end embed in the thermostat 
with temperature T_ = 0.9T. In this case the dynam- 
ics of the atoms with numbers n = 1, ...,Ne is described 
by Langevin equations (46) with temperature T = T+, 
the atoms with numbers n — N — + 1, N are de- 
scribed by these equations with temperature T = T_, 
while the atoms in the central part of the nanoribbon, 
with numbers n = + 1, N — N^, are described by 
the frictionless equations, Eqs. (46) with F = and 
rjkn = 0. We take the relaxation time tr = 100, at which 
the edge effects at the interfaces between the ribbon and 
the thermostats are not essential. We integrate the equa- 
tions of motion to find the stationary energy fiux J along 
the ribbon. 

To determine the temperature profile, first we find en- 
ergy per mode at the longitudinal coordinate x = na hy 
averaging the distribution across the ribbon of the scalar- 
model double kinetic energy: 



K N K N-l K 



k=l n=l 



k=l n=l 



fc=l 



(45) 

where the dot denotes the derivative with respect to the 
dimensionless time t. 

For the dimensionless quantities, equations of motion 
will take the form: 



Ukn 



an 

dUkn 



-Tiikn+Vkn, k^l,...,K, n = l,...,iV 

(46) 

where F = Tto is the dimensionless friction, and the di- 
mensionless random forces rjkn is the color noise with the 



1 



Bki 



K 



(49) 



Since in the ribbon with a constant average width the 
phonon spectrum does not depend on the longitudinal 
coordinate and the inter-mode interaction is absent in 
the harmonic system, we can determine the local temper- 
ature from the equation ekin{n) = e{Tn)^ similar to Eq. 
(39). Figure 8 shows the distribution along the ribbon of 
energy fiux J, kinetic energy per mode Ckin and temper- 
ature T. There is a stationary flux in the area between 
the thermostats. Since the temperature distribution of 
phonon frequencies in the ribbon with harmonic inter- 
actions is determined by the Langevin thermostats with 
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FIG. 9: (Color online) (a) Thermal conductivity k of rough- 
edge ribbon (ribbon width K = 12, rough edges widths K\ = 
4, porosity of rough edges p = 0.05) versus temperature T for 
ribbon length iV = 600 (line 1) and iV = 400 (line 2). Dashed 
lines 3 and 4 give the dependencies for perfect nanoribbons 
(with zero porosity of edges, p = 0) with A*" = 600 and = 
400, respectively. Inset shows the low-temperature limit, T < 
0.22. (b) Factor of thermal conductivity suppression by rough 
edges r versus temperature for ribbon length A'' = 600 (line 
5) and N = 400 (line 6). 



Bose-Einstein color noise, we can uniquely determine the 
local temperature T„, which monotonously changes from 
the value at the left end to the value r_ at the right 
end with a linear gradient, see Fig. 8(c). 

Then the dimensionless thermal conductivity of the 
finite-length ribbon can be found as 

ti{N,T) ^{N~ 2Ne)J/[K{T+ - T_)]. (50) 

In the ribbon with perfect (atomically smooth) edges, 
phonons have infinite mean free path and therefore the 
energy flux J does not depend on the length of the cen- 
tral part of the ribbon Nc = N — 2N^. In this case, as 
one can conclude from Eq. (50), the thermal conductiv- 
ity increases linearly with Nc and therefore such ribbon 
behaves as an ideal (ballistic) thermal conductor. In the 
absence of anharmonicity of interatomic interactions, this 
conclusion is valid for any temperature. 

We consider a ribbon which consists oi K = 12 
molecular chains. To model the roughness of the rib- 
bon edges, we delete with probability p = 1 — d some 
atoms from the chains with numbers k = l,...,Ki and 
k = K — Ki + 1, K . Here Ki is a width of the rough 
edges, < d < 1 is their atomic density and the parame- 
ter p determines in fact the porosity of the ribbon lattice 



in the defect edges. We take in the following Ki = 4 
and p = 0.05. We computed the thermal conductiv- 
ity k{N,T) for N = 400 {Nc = 200) and = 600 
{Nc = 400), see Fig. 9. As one can see from this fig- 
ure, at high (dimensionless) temperature the ribbon with 
the length N = 600 has lower thermal conductivity than 
the ribbon with length TV = 400. This means that the 
ribbon behaves much like thermal insulator, in which the 
thermal conductivity decreases with the length of the sys- 
tem, see also [11]. But the situation is changed for low 
temperature, for T < 0.18, when the longer ribbon has 
the higher thermal conductivity, as in the case of perfect 
nanoribbons, see the inset in Fig. 9 (a). 

In this connection we can define the factor of ther- 
mal conductivity suppression by rough edges r{T, N) = 
Ki{N, T)/ko{N, T), where ko{N, T) is thermal conductiv- 
ity of a ribbon with ideal (atomically smooth) edges and 
length N , ki{N, T) is thermal conductivity of the rough- 
edge ribbon with the same length and width. As one can 
see from Fig. 8 (b), the rough edges suppress the ther- 
mal conductivity for all temperatures, r < 1 for T > 0, 
but the suppression monotonously decreases with the de- 
crease of temperature: for T — >• the thermal conductiv- 
ities of the ideal- and rough-edge ribbons flatten, r — >• 1. 
This means that long wave acoustic phonons are not scat- 
tered by surface roughness and therefore at low enough 
temperature the rough-edge quantum ribbon becomes an 
ideal (ballistic) thermal conductor, in which the thermal 
conductivity increases with the conductor length. This 
in turn means that at low enough temperature we ap- 
proach the limit of ballistic quantum thermal transport, 
when the value of thermal conductance Gth = J / ^T, 
the ratio of the thermal flux J through the quasi-one- 
dimensional thermal conductor to the temperature dif- 
ference AT = T+ —T-, has a quantized value TT^kgT/3h 
(per each of the four massless acoustic modes of the quasi- 
one-dimensional waveguide), which does not depend on 
the material properties (and perfectness) of the thermal 
conductor, see, e.g., Refs. [17-19]. This property of low- 
temperature thermal transport is in drastic contrast to 
that in the classical high-temperature regime, in which 
the same quasi-one-dimensional rough-edge nanoribbon 
behaves as a thermal insulator, in which the thermal 
conductivity decreases with the insulator length, sec Ref. 
[11]. 



VII. MODELING OF HEAT TRANSPORT IN 
CARBON NANOTUBE 

The computation of the thermal conductivity can be 
performed by two methods. In the first method, the sys- 
tem is thermalized with the use of equilibrium molecular 
dynamics and then the interaction with the thermostat 
is switched off. The current-current correlation function 
in the system is computed and the coefficient of ther- 
mal conductivity is found with the help of Green-Kubo 
formula based on this correlation function. In the sec- 
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ond method, the method of non-equihbrium dynamics, 
the direct modehng of thermal transport is performed. 
For this purpose, two ends of the quasi-one-dimensional 
system are embedded in two thermostats with different 
temperatures. The stationary flux of energy is computed 
and the coefficient of thermal conductivity is determined 
from the known energy flux, temperature difference and 
length of the system. In the approach based on the non- 
equilibrium dynamics, the correct use of thermostats is 
important. The Langevin thermostat can be used in such 
approach. 

The both methods model the heat transport within 
the classical Newtonian dynamics and therefore do not 
allow to take into account the quantum effects. These 
methods can be justified only for high enough tempera- 
ture when the statistics of the system becomes completely 
classical. Nevertheless the method of the non-equilibrium 
dynamics can be easily adapted for the use in the semi- 
quantum approach to molecular dynamics modeling of 
thermal transport. Below we demonstrate this by exam- 
ple of a carbon nanotube. 

For the classical modeling of heat transport, we con- 
sider finite nanotube with chirality index (5,5) (see 
Fig. 4), and embed its left end in Langevin thermostat 
with temperature T+ = I.IT and the right end in a ther- 
mostat with temperature T_ = 0.9T, where T is a tem- 
perature at which the thermal conductivity will be de- 
termined. For this purpose the dynamics of the atoms at 
the left (right) end of the nanotube should be described 
by the Langevin equations (1) with a white noise with the 
correlation function given by Eq. (5) with temperature 
T = T_|_ (T_). Dynamics of the atoms in the central part 
of the nanotube we describe with the Hamilton equations 
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FIG. 10: (Color online) Distribution of (a) local energy flux J 
and (b) local temperature T along the carbon nanotube (5,5) 
C6020, X axis is along the nanotube axis. Nanotube length is 
L = 75.04 nm. Gray arias indicate the nanotube ends, em- 
bedded in classical Langevin thermostats with temperatures 
r± = (1 ± (5)300K, for 5 = 0; 0.0125; 0.025; 0.05; 0.1 (curves 
1 and 6; 2 and 7; 3 and 8; 4 and 9; 5 and 10, respectively). 



-dH/dr„ 



(51) 



This allows us to find the distribution of temperature 
T{x) and energy flux J{x) along the nanotube, see 
Fig. 10. The detailed calculation of these quantities can 
be found in Ref. [13]. In the classical molecular dynamics 
description, temperature of a particle can be determined 
through its average double kinetic energy per mode (or 
per degree of freedom) as 



T = M{X^ +f + Z^)/3kB = CMn/ki 



(52) 



where M and {x^y, z) are particle mass and cartesian 
coordinates. 

To illustrate the modeling of heat transport, we con- 
sider a nanotube with a flxed length L = 75.04 nm (which 
corresponds to N=6020 atoms) and embed its ends with 
lengths Le = L/3 in classical Langevin thermostats with 
temperature T± = (1 ± 5)300 K, where 6 = 0, 0.0125, 
0.025, 0.05, 0.1. For the relaxation time, we take tr = 0.4 
ps. As one can see in Fig. 10, for S = when the 
ends have the same temperature = T- = 300 K, 
the temperature is a constant along the whole nanotube 
{T{x) = 300 K) and there is no heat flux in the circuit 
{j{x) = 0). For T+ > T_ (5 > 0), the linear temperature 



gradient T{x) and heat flux are formed in the central part 
of the nanotube: J{x) = J > for < x < L — L^- 
This allows to determine the thermal conductivity of the 
nanotube of length L as 



k(L,T) = (L- 2Le)J/5'(T+ -T_), 



(53) 



where S 



is nanotube cross section [radius r is 



equal to 0.335 nm for the (5,5) single- walled nanotube]. 
One can obtain more accurate estimate for the thermal 
conductivity from the temperature profile in the central 
part of the nanotube: 

-K={L- 2Le)J/S[T{Le) - T{L - L^)], (54) 

where T(x) is a distribution of the particle kinetic en- 
ergy (temperature) along the nanotube. The definition 
of thermal conductivity, given by Eq. (54), allows one 
to separate boundary resistances from the thermal resis- 
tance of the nanotube by itself. The boundary (Kapitza) 
resistance causes finite difference between the tempera- 
tures in the bulk of the heat reservoir and just at the 
interface with the nanotube, see Fig. 10(b). Since one 
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has T{Le) < T+ and T{L - L^) > T_ at the inter- 
faces between the ends embedded in the thermostats and 
the central part of the nanotube, Eqs. (54) always give 
higher values of thermal conductivity than Eq. (53) does: 
K > K. Dependencies of the obtained values of k and R on 
the temperature difference at the ends of the nanotube, 
AT = 25 X 300 K, are presented in Table II. 

As one can see in Table II, for AT < 60 K the halving 
of the temperature difference results in the halving of the 
energy flux J. Since in this case the flux is proportional 
to the temperature difference, the determination of the 
thermal conductivity k (k) does not depend on the value 
of AT. It is worth mentioning that the accuracy of the 
determination of the energy flux J decreases with the 
decrease of AT. Therefore the most convenient for the 
MD modeling of thermal conductivity are the values of 
the temperatures at the nanotube ends as T+ = I.IT 
and T_ = 0.9T, although a relatively high gradient of 
temperature is formed in this case (0.8 K/nm for T = 
300 K). 

For the semi-quantum MD modeling of heat trans- 
port, we must use the Langevin equations (1) with 
the random forces ^q,„ with the power spectral den- 
sity, given by {^an£,pm)ui = 2MTkBT{SianSi 
Eq. (34), with temperature T = r+ or T = T_. The 
color noise is determined as ^Q„(t) — Can{^T/kBT±) = 
kBT±^y 2MT I hSian {t) . The zero-point oscillations do 
not contribute to the thermal transport, and therefore 
they will not be taken into account in the following and 
in all the formulas below we will put S'oan(''") = 0. The 
thermal conductivity k can be found with the use of Eq. 
(53). 

In the semi-quantum MD modeling, the mean value of 
particle kinetic energy ekin does not coincide with the 
temperature. This value depends not only on the tem- 
perature, but also on the density of plionon states. For 
example, the carbon atoms at the end hemispheres of 
the nanotube have different vibrational spectrum than 
the atoms in the central part of the nanotube. As one 
can see in Fig. 11 (b), at the very ends of the nanotube, 
for X close to and to 75 nm, the mean kinetic energy of 
carbon atoms is higher than that in the central part of 
the nanotube. In the central part of the nanotube, all the 
atoms have the same vibrational spectrum. Here Cfcm de- 
pends unambiguously on temperature and therefore the 
distribution of kinetic energy along the nanotube charac- 
terizes unambiguously the distribution of temperature. 

Carbon nanotube is a stiff molecular system with 
weakly nonlinear dynamics. One can compute its vibra- 
tional spectrum and then use Eq. (39) to determine the 
local temperature in the system. 

We consider the nanotube of length L = 75.04 nm, 
which ends with the length Le = L/3 are embedded 
into semi-quantum Langevin thermostats with temper- 
atures T± = (1 ± 5)300 K (when 6 = 0, 0.0125, 0.025, 
0.05, 0.1). As one can see in Fig. 11 (b), for r+ = 
T_ = 300 K {6 = 0) the mean particle kinetic en- 
ergy Ckin is not constant along the nanotube. In the 
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FIG. 11: (Color online) Distribution of (a) local energy flux 
J{x), (b) local kinetic energy per mode ekin{x), and (c) lo- 
cal temperature T{x) along the carbon nanotube (5,5) C6020 
axis X. Nanotube length is L = 75.04 nm. Gray arias indi- 
cate the nanotube ends, embedded in semi-quantum Langevin 
thermostats with temperatures T± = (1 ± 5)300K, for 5 = 0; 
0.0125; 0.025; 0.05; 0.1, curves 1 and 6; 2 and 7; 3 and 8; 
4 and 9; 5 and 10, respectively. Thin lines show distribu- 
tions of average particle energies eki„{x) in the central part 
of the nanotube without the artefact surplus thermal energy 
Aefei„(x). Local temperature T{x) in (c) is determined from 
Eq. (39) with the use of ekin{x), given by thin lines in (b). 



central part, which does not interact directly with the 
thermostats, particles have higher thermal energy than 
in the end parts, which do interact directly with the 
thermostats. For the end parts the mean kinetic en- 
ergy is ekin{T±) = 53.5 K, while in the central part 
one has ekin{T±) < ekin{x) < 55.5 K. The surplus ki- 
netic energy Aefc„(a;) = ekin{x) - efci„(T±), which is 
< Aekinix) < 2K, is an artefact of the semi-quantum 
molecular dynamics approach when the latter is applied 
to nonlinear lattice, but it does not produce any thermal 
flux: J{x) = for Le < X < L — Le, see Fig. 11 (a). 
The appearance of the surplus kinetic energy in the 
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TABLE 11: Heat flux J, thermal conductivities k and k of the (5,5) carbon nanotube of length L — 75.04 nm versus temperature 
difference at the ends AT — T+ — T- [mean temperature is T = (T+ +T_)/2 = 300 K, lengths of the embedded in thermostats 
ends are Le — L/3, relaxation time is tr — 0.4 ps]. The values were obtained with the use of classical molecular dynamics 
(CMD) and of semi-quantum molecular dynamics (SQMD). The values of k and R have been calculated with the use of Eqs. 
(53) and (54). 





CMD 


SQMD 


AT (K) 


J (eVA/ps) 


K (k) (W/mK) 


J (eVA/ps) 


K (k) (W/mK) 


60.0 


1.714 


260.0 (448) 


0.659 


99.9 (260) 


30.0 


0.856 


259.5 (447) 


0.315 


95.4 (249) 


15.0 


0.426 


258.6 (443) 


0.154 


93.4 (244) 


7.5 


0.210 


257.1 (436) 


0.077 


93.0 (243) 



nonlinear system is related with phonon interaction 
caused by the anharmonicity, which induces the trans- 
fer of energy from the low-frequency modes to the high- 
frequency ones. In the nanotube ends, such transfer is 
suppressed by viscous friction, but in the central part of 
the nanotube it can result in the surplus excitation of the 
high-frequency modes. In the case of the white-noise heat 
baths, such inter-mode transfer does not occur because of 
the equipartition excitation. The inter-mode transfer re- 
sults in surplus excitation of high-frequency modes in the 
central part of the nanotube, which ends are embedded in 
color-noise thermostats. The analysis of the vibrational 
spectrum shows that weakly anharmonic interparticle po- 
tential in the carbon nanotube results in a slight accu- 
mulation of energy in the high-frequency modes, with 
uj > 900 cm~^, while the low-frequency modes, with 
uj < 900 cm^^, follow the Bose-Einstein distribution, see 
Fig. 12. As it was explained above in connection with the 
quantum definition of the lattice temperature, see Eqs. 
(14) and (17), this means that the high-frequency modes 
possess higher effective temperature (and therefore are 
"hot phonons") then the lattice temperature, which is 
determined in turn by the Bose-Einstein distribution of 
the low- frequency modes. 

For the system with the harmonic interparticle poten- 
tial, the Bose-Einstein distribution is valid for all phonon 
modes in all the lattice system with the ends, embedded 
in color-noise thermostats. Indeed, as we have shown in 
Section V, there is no surplus excitation of high-frequency 
modes in the central part of the nanoribbon with the 
harmonic interparticle potential. 

For > T- ((5 > 0), a constant heat flux is formed 
in the central part of the nanotube {J{x) = J > for 
Lg < X < L ~ Le), which is proportional to the tem- 
perature difference AT = — see Table II. As 
one can see in Fig. 11 (b), the distribution of local ki- 
netic energy per mode ekin{x) has a nonlinear form for 
Le < X < L—Lp - But after the subtraction of the artefact 
surplus contribution Ackinix) from the average kinetic 
energy ekin{x), the remaining local energy Ckinix) has 
a linear slope. If one determines the local temperature 
T(x) from the distribution of the local energy ekin{x) 
with the use of Eq. (39), one gets the linear distribu- 
tion of local temperature along the tube, see Fig. 11 (c). 



The same linear distribution of local temperature can 
also be obtained from Eq. (40), in which the integra- 
tion over frequencies is performed in the low-frequency 
domain wi = < cj < W2 = 900 cm^^. As one can 
see in Fig. 12, in this frequency domain the vibrational 
spectrum has a correct Bose-Einstein distribution, which 
allows the correct definition of temperature with the use 
of Eq. (39) in all parts, embedded ends and free central 
part, of the nanosystem. The linear distribution of local 
temperature along the tube allows in turn to determine 
more accurately the thermal conductivity of the system 
with the use of Eq. (54). 

As one can see from Table II, for AT < 60 K the value 
of thermal conductivity changes only weakly with the 
change of temperature difference. Thus for the numerical 
modeling, one can use T± = (1 ± 8)T with 5 = 0.1. Such 
temperature difference allows one to find rather fast the 
distribution of the heat flux and temperature along the 
nanotube, and the obtained value of the thermal conduc- 
tivity K differs little from the values obtained for smaller 
5. 

Now we turn to the analysis of the vibrational fre- 
quency spectrum of carbon atoms near the nanotube 
ends embedded in thermostats with different tempera- 
tures, r+ = 330 K and r_ = 270 K, and in the nan- 
otube central part. As one can see in Fig. 12, the power 
spectral density of thermal vibrations at the nanotube 
ends corresponds to the quantum statistics of phonons. 
But there are some distinctions in the nanotube central 
part, namely the high-frequency modes with frequencies 
ijj > 900 cm^^ are excited more strongly. This effect 
is related with the phonon interaction caused by anhar- 
monicity, which induces the transfer of energy from the 
low-frequency vibrations to the high-frequency ones. In 
the nanotube ends, such transfer is suppressed by vis- 
cous friction, but in the central part of the nanotube it 
can result in the surplus (artefact) excitation of the high- 
frequency modes. As one can see from Fig. 12 (b), this 
effect is rather weak. It does not change the heat flux 
considerably for the used values of the nanotube length. 
It is worth mentioning in this connection that the high- 
frequency modes, which are revealed in Fig. 12 (b) for 
Lo > 900 cm~^, determine, via Eq. (14) or Eq. (17), 
the effective temperature of hot phonons, which is higher 
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FIG. 12: (Color online) Power spectral density of thermal 
vibrations in carbon nanotube (5,5) with length L — 75.04 nm 
(a) for the atoms at the left end with temperature 7+ — 
330 K (curve 2, red), right end with temperature T- = 270 K 
(curve 1, blue), (b) Line 3 (red) gives vibrational spectral 
density in the central part of the nanotube, which does not 
interact with the thermostats, line 4 (blue) gives the spectral 
density given by Eq. (37) for the corresponding system of 
quantum oscillators. The frequency spectrum of the mean 
kinetic energy per atom is shown. 

than the temperature of the equilibrium, Bose-Einstein, 
phonons. But the surplus hot-phonon temperature, sim- 
ilar to the surplus kinetic energy of the particles in the 
central part of the nanotube, shown in Fig. 11 (b), does 
not contribute to the linear distribution of the temper- 
ature and therefore to the thermal conductivity of the 
nanotube. 

The difference in the power spectral density close to 
zero frequency in the end, Fig. 12 (a), and central. Fig. 
12 (b), parts of the nanotube is related with the use of 
the effective fix-end and free-end boundary conditions in 
the simulation of corresponding spectra. 

We would like to mention that the obtained results for 
the temperature dependence of the specific heat and ther- 
mal flux in the considered carbon nanotubes almost do 
not depend on the considered relaxation times = 0.1, 
0.2, 0.4, 0.8 ps. The longer relaxation time increases the 
computation time. The steady-state heat transport es- 
tablished longer time and the computation of the mean 
values requires the integration along longer phase-space 



trajectories. The shorter relaxation time the effective 
viscosity smears out vibrational spectra of particles in- 
teracting with the thermostats. It can also increase the 
reflectivity of high-frequency phonons at the interfaces 
between the nanotube central part and the thermostats. 
The optimal relaxation time for the carbon nanotube is 
U = 0.4 ps. 

Now we consider the nanotube of the fixed length 
L = 50.08 nm (which corresponds to = 4020 atoms), 
with Le = i/4 and t,. = 0.4 ps. Temperature de- 
pendence of the thermal conductivity of such nanotube 
is shown in Fig. 13. As one can see from this figure, 
within the classical description the thermal conductivity 
Ki monotonously increases with the decrease of temper- 
ature. This property of "classical" thermal conductiv- 
ity is related with the decrease of anharmonicity of lat- 
tice dynamics with the decrease of temperature, which 
in turn results in the increase of phonon mean free path. 
But the situation changes drastically with an account 
for quantum statistics of phonons. In the semi-quantum 
description, thermal conductivity K2 first increases with 
the decrease of temperature from the high enough one 
but then it reaches its maximal value at T = 350K and 
monotonously decreases to zero (k2 for T — 0). 
Such temperature dependence of thermal conductivity, 
which is characteristic for solids [15, 16], is related with 
the quantum decrease of the specific heat of the system. 
This is confirmed by the property that at low temper- 
ature T < 250K the thermal conductivity of the nan- 
otube in the semi-quantum approach K2{T) is described 
with high accuracy by the relation k.2{T) ~ c{T)ki{T), 
where ki{T) is thermal conductivity of the carbon nan- 
otube computed within the classical description, c{T) = 
C{T)/3NkB is nanotube dimensionless specific heat at 
temperature T, see Figs. 5 and 13. 

Length dependence of the nanotube thermal conduc- 
tivity is shown in Fig. 14. As one can see from this fig- 
ure, both methods give monotonous increase of the ther- 
mal conductivity with the nanotube length. The semi- 
quantum description gives the lower value of the conduc- 
tivity than the classical one for all the considered nan- 
otube lengths. As one can see from this figure, line 2 
lies below line 3 for relatively short nanotubes, which is 
consistent with Fig. 13. But at some nanotube length 
the line 2 intersects the line 3. This means that for 
longer nanotube lengths the mean free path of "semi- 
quantum" phonons is larger than that of purely " classi- 
cal" phonons. This in turn is related with the decrease 
in general of phonon mean free path with phonon fre- 
quency and the property that the mean frequency of semi- 
quantum phonons is always lower than that of classical 
phonons. The intersection at some nanotube length of 
the line 2 with line 3 in Figs. 13 and 14 means that short 
enough nanotubes effectively filter out the high-frequency 
and short-mean-free-path phonons, as it occurs in meso- 
scopic one-dimensional samples, see, e.g., Ref. [17]. Line 
2 should intersect the line 1 at even longer lengths, when 
thermal conductivity of the "semi-quantum" nanotube 
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FIG. 13: (Color online) Temperature dependence of thermal 
conductivity Ki, (i = 1,2), of the carbon nanotube (5,5) 
with length L — 50.08 nm: solid lines 1 and 2 are ob- 
tained within the classical and semi-quantum descriptions, 
respectively; dashed line 3 gives the temperature dependence 
c(r)«:i(r), where c{T) = C{T)/3NkB is dimensionless spe- 
cific heat of the carbon nanotube at temperature T. 

will become larger than that of purely "classical" nan- 
otube. Figure 15 (b) below shows that similar situation 
can also be realized in a nanoribbon. Concerning carbon 
nanotubes, atT = 300K such intersection can be realized 
only in very long nanotubes because of their very high 
Debye temperature (which in fact exceeds the melting 
temperature of the material). 

VIII. THERMAL CONDUCTIVITY OF 
NANORIBBON WITH PERIODIC 
INTERATOMIC POTENTIALS 

In order to demonstrate the characteristic tempera- 
ture dependence of thermal conductivity of quantum low- 
dimensional system with strongly nonlinear interatomic 
interactions, we consider a nanoribbon with periodic 
interatomic potentials and perfect (atomically smooth) 
edges, cf. Eqs. (41)-(44): 

V{r) ei[l - cos(r)], U{r) = 62(1 - cos(r)] . (55) 

In the following we will consider ei = 1, £2 = 0.5, when 
for small relative displacements of the nearest-neighbor 
atoms r the potentials (55) will coincide with the har- 
monic potentials given by Eq. (44). With the periodic 
interatomic potentials, given by Eq. (55), the ribbon be- 
comes a system of coupled rotators. 

A chain of coupled rotators presents a unique 
translationally-invariant one-dimensional system with a 
finite thermal conductivity [20, 21]. In such system roto- 
breathers can be excited, which cause strong momentum- 
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FIG. 14: (Color online) Thermal conductivity k versus length 
of the central part Lc = L — 2Le of the carbon nanotube 
(5,5) obtained with the use of classical (solid line 1) and semi- 
quantum (solid line 2) descriptions, and thermal conductivity 
c(r)fi:i(r) (dashed line 3, as line 3 in Fig. 11) at T = 300K. 
The length Lc is measured in A, the thermal conductivity k 
- in W/mK. 



nonconserving scattering of phonons and result in finite 
thermal conductivity of the translational-invariant sys- 
tem. The density of rotobreathers increases with temper- 
ature increase and correspondingly the phonon mean free 
path decreases. In the classical description, this causes 
the monotonous decrease of thermal conductivity k with 
the increase of temperature: k — > for T — ^ 00. In 
the opposite limit T — >■ 0, the phonon mean free path 
monotonously increases and therefore the thermal con- 
ductivity diverges: k — >• 00 for T — >■ 0. 

We consider a ribbon with a width K — 2 and length 
N = 800. We embed the ribbon ends with length = 
300 each into Langevin thermostats with the color noise 
with temperature T+ = I.IT and T- = 0.9T in the left 
and right edge, respectively. For the relaxation time, we 
take tr = 100. Integration of Eqs. (46) for n = 1, ...,iVe 
and n^N — Ne + l,.--,N, and of frictionless equations, 
Eqs. (46) with f = and ijkn, for n = Ne + 1, ...,N-Ne, 
allows us to find an average energy flux along the ribbon 
J. With the use of Eq. (50), we then obtain the thermal 
conductivity of the finite-length nanoribbon. 

Results of numerical modeling are presented in Fig. 15. 
As one can see from the temperature dependence of the 
nanoribbon specific heat. Fig. 15 (a), within the semi- 
quantum description the anharmonicity of atomic dy- 
namics starts to show up for T > 0.5. Just at this 
temperature specific heat of the ribbon with the non- 
linear interatomic potentials (55) starts to deviate from 
the specific heat of the ribbon with harmonic potentials 
(44), see Fig. 15 (a), line 2. It is worth mentioning that 
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FIG. 15: (Color online) (a) Specific heat c of an ideal-edge 
ribbon with periodic interatomic potential (55) versus tem- 
perature T, obtained within the classical (line 1) and semi- 
quantum (line 2) descriptions. Dashed lines 3 and 4 show 
the dependencies for the nanoribbon with the harmonic in- 
teratomic potential (44) within, respectively, the classical and 
semi-quantum descriptions, (b) Thermal conductivity k. ver- 
sus temperature T for the ribbon with length A'^ = 800 and 
width K = 2 (the length of each of the two ends embedded 
into the thermostats with different temperatures is A^e ~ 300): 
lines 5 and 6 give the results obtained within the classical and 
semi-quantum descriptions, respectively. 



in the classical description the anharmonicity of lattice 
dynamics starts to show up for lower temperatures than 
that in the semi-quantum approach. This is related with 
the property that the anharmonicity is more pronounced 
in the dynamics of short wave phonons and therefore the 
quantum freezing out of the high frequency oscillations 
results in the effective decreasing of the dynamics anhar- 
monicity (nonlinearity) at low temperature. 

In the classical description, the thermal conductivity 
of the ribbon with nonlinear interatomic potentials (55) 
decreases monotonously with the increase in tempera- 
ture: K — !• for r — > oo, see Fig. 15 (b), line 5. Within 
the semi-quantum description, the ribbon thermal con- 
ductivity reaches its maximal value at T = 0.6. For the 
lower and higher temperatures, the thermal conductivity 
decreases, k — > for T — and k ^> for T — oo, see 



Fig. 15 (b), line 6. For the low temperature T < 0.5, 
system dynamics becomes almost linear when phonons 
propagate along the ribbon almost ballistically and the 
decrease of thermal conductivity is related with the quan- 
tum decrease of the specific heat. In this low temperature 
limit the thermal conductance of a short enough ribbon, 
in which phonon mean free path is longer than the ribbon 
length, can reach the lowest (quantum) value which has 
a linear temperature dependence, see line 6 in Fig. 15 (b) 
and compare the corresponding limit for the rough-edge 
ribbons, shown in Fig. 9, (a) and (b). 

For temperature T > 0.4, the semi-quantum descrip- 
tion gives higher value for the thermal conductivity than 
that of the classical description. This is also related with 
the quantum freezing out of the high-frequency oscilla- 
tions, which results in the decrease of the mean frequency 
of thermal phonons and correspondingly in the increase 
of phonon mean free path and phonon thermal conduc- 
tivity. Therefore for T > 0.4 the effect of the increase 
of phonon mean free path exceeds the effect of the de- 
crease of low temperature specific heat. At high tem- 
perature, when phonons can be described with the clas- 
sical statistics, the thermal conductivities given by the 
classical and semi-quantum descriptions merge, but the 
conductivity given by the semi-quantum description is al- 
ways higher because of higher phonon mean free path un- 
der the equal, classical and semi-quantum, specific heats. 
The temperature dependence of thermal conductivity, 
given by line 6 in Fig. 15 (b), reminds the known temper- 
ature dependence of thermal conductivity in insulators, 
when the temperature of maximal thermal conductiv- 
ity separates the low-temperature boundary-scattering 
regime from the high-temperature anharmonic Umklapp- 
scattering regime, see, e.g., Refs. [15, 16, 22, 23]. In the 
case of the nanoribbon with periodic interatomic poten- 
tials and perfect edges, the maximal thermal conductiv- 
ity, clearly displayed by curve 6 in Fig. 15 (b), is reached 
for the temperature, at which phonon mean free pass, 
which is determined in this system by the anharmonic 
scattering, levels off with the ribbon length. For the 
higher temperature, the phonon mean free pass becomes 
shorter then the ribbon length. Figure 15 (b) presents 
one of the main results of this work. 



IX. SUMMARY 

In summary, we present a detailed description of semi- 
quantum molecular dynamics approach in which the dy- 
namics of the system is described with the use of clas- 
sical Newtonian equations of motion while the effects of 
phonon quantum statistics arc introduced through ran- 
dom Langevin-like forces with a specific power spectral 
density (the color noise). We describe the determination 
of temperature in quantum lattice systems, to which the 
equipartition limit is not applied. We show that one can 
determine the temperature of such system from the mea- 
sured power spectrum and temperature- and relaxation- 
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rate-independent density of vibrational (phonon) states. 
We have applied the semi-quantum molecular dynam- 
ics approach to the modeling of thermal properties and 
heat transport in different low-dimensional nanostruc- 
tures. We have simulated specific heat and heat trans- 
port in carbon nanotubes, as well as the heat trans- 
port in molecular nanoribbons with perfect (atomically 
smooth) and rough (porous) edges, and in nanoribbons 
with strongly anharmonic periodic interatomic poten- 
tials. We have shown that the effects of quantum statis- 
tics of phonons are essential for the carbon nanotube 
in the whole temperature range T < 500K, in which 
the values of the specific heat and thermal conductiv- 
ity of the nanotube are considerably less than that ob- 
tained within the description based on classical statis- 
tics of phonons. This conclusion is also applicable to 
other carbon-based materials and systems with high De- 
bye temperature like graphene, graphene nanoribbons, 
fullerene, diamond, diamond nanowires etc. We have 
shown that quantum statistics of phonons and porosity 
of edge layers dramatically change low temperature ther- 
mal conductivity of molecular nanoribbons in compar- 
ison with that of nanoribbons with perfect edges and 
classical phonon dynamics and statistics. The semi- 
quantum molecular dynamics approach has allowed us to 
model the transition in the rough-edge nanoribbons from 
the thermal-insulator-like behavior at high temperature, 
when the thermal conductivity decreases with the con- 
ductor length, see Ref. [11], to the ballistic-conductor- 
likc behavior at low temperature, when the thermal con- 
ductivity increases with the conductor length. We have 
also shown that the combination of strong nonlinearity 



of the interatomic potentials with quantum statistics of 
phonons changes drastically the low-temperature ther- 
mal conductivity of the system. The thermal conductiv- 
ity in such samples demonstrates very pronounced non- 
monotonous temperature dependence, when the temper- 
ature of maximal thermal conductivity separates the low- 
temperature ballistic phonon conductivity from the high- 
temperature anharmonic-scattering one. At the temper- 
ature of maximal thermal conductivity, phonon mean 
free pass levels off with the length of the perfect-edge 
anharmonic quasi-one-dimensional system. Such non- 
monotonous temperature dependence of thermal conduc- 
tivity is known in bulk insulators and is very different 
from monotonously decreasing with temperature conduc- 
tivity of nanoribbons with the same interatomic poten- 
tials and classical phonon dynamics and statistics, cf. 
Refs. [20, 21]. 
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